МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ
НУ ”Львівська політехніка”
Звіт
про виконання лабораторної роботи №3
з курсу: ”Алгоритмічні основи криптології”
на тему:
”Метод Ньютона для розв’язування систем нелінійних рівнянь ”
Мета роботи: ознайомлення з найпоширенішим ітераційним методом розвязування систем нелінійних рівнянь- методом Ньютона.
Теоретичні відомості
Систему нелінійних рівнянь у загальному випадку можна зобразити у вигляді
(1)
тобто як n функцій від n невідомих причому функції не обов’язково лінійно залежать від зміних .
Позначимо вектор зміних через , а вектор функції через . Тоді систему можна запиати у формі
Модифікований метод Ньютона
При використані стандартного методу Ньютона на кожній ітерації доводиться обчислювати новий якобіан , хоч зрозуміло, що при закінчені ітерації він повинен прийняти стабільне значення , - розвязок. У модифікованому або спрощеному або спрощеному методі Ньютона якобіан заміняють правильно підібраною матрицею А. Звичайно, найкращим, але практично недосяжним була б заміна , - розвязок.
Але на практиці користуються компромісним рішенням:
-вибирають за якобіан в початковій точці , а ітерації проводяться за наступною формулою
;
-зберігають А протягом певного числа ітерацій;
- на певній r-й ітерації змінюють А, прирівнюючи її якобіану і з новим значенням знову виконують певне число ітерацій і т.д.
Отже, якобіан обчислюється тільки час від часу, за рахунок чого досягається економія машиного часу. Однак, збіжність методу при цьому стає практично лінійною.
Завдання: Розв’язати систему нелінійних рівнянь модифікованим методом Ньютона, вибираючи за початкові наближення . Ітерації проводити до збігу двох послідовних наближень з похибкою . Обертання матриці здійснити методом LU-факторизації, похідні методом кінцевих різниць.
Текст програми
#include <stdio.h>
#include <conio.h>
#include <math.h>
#define n 2
#define e 0.0001
#define h 10e-8
float func(float *a,int i)
{
if (i==1) return (a[1]-cos(a[2])-1.0);
else return (a[2]-sin(a[2])-1.0);
}
main()
{
float x[n+1],x1[n+1];
float max, xp,yp;
int i,k,f,j,m;
float u[n+1][n+1],l[n+1][n+1],Jcob[n+1][n+1];
float Jc[n+1][n+1];
float s1=0, s2=0;
float b[n+1],y[n+1];
x[1]=x1[1]=01.0;
x[2]=x1[2]=01.0;
do
{
xp=x[1];
yp=x[2];
/*----Metod kincevyh riznyts---*/
for (i=1;i<=n;i++)
for (j=1;j<=n;j++)
{
x1[j]+=h;
Jc[i][j]=(func(x1,i)-func(x,i))/h;
x1[j]-=h;
}
/*----------LUrozklad i obern------------------------------*/
for (f=1;f<=n;f++)
{
for (i=1;i<=n;i++)
{
if (f==i) b[i]=1;
else b[i]=0;
}
/*Vyznachennia elementiv L ta U*/
for (i=1;i<=n;i++)
l[i][1]=Jc[i][1];
for (j=2;j<=n;j++)
u[1][j]=Jc[1][j]/Jc[1][1];
k=2;
i=k;
for(m=1;m<=(k-1);m++)
s1+=(l[i][m]*u[m][k]);
l[i][k]=Jc[i][k]-s1;
s1=0;
y[1]=b[1]/l[1][1];
for(i=2;i<=n;i++)
{
for(m=1;m<=(i-1);m++)
s2+=(l[i][m]*y[m]);
y[i]=(b[i]-s2)/l[i][i];
s2=0;
}
x1[n]=y[n];
for(i=n-1;i>=1;i--)
{
for (m=i+1;m<=n;m++)
s2+=u[i][m]*x1[m];
x1[i]=y[i]-s2;
s2=0;
}
//------Pereprysvoennia--
if (f==1)
{
Jcob[1][1]=x1[1];
Jcob[2][1]=x1[2];
}
else
{
Jcob[1][2]=x1[1];
Jcob[2][2]=x1[2];
}
}
for (i=1;i<=n;i++)
{
s1=0;
for (j=1;j<=n;j++)
{
Jcob[i][j]*=func(x,j);
s1+=Jcob[i][j];
}
if (i==1) x[i]=x1[i]=xp-s1;
else x[i]=x1[i]=yp-s1;
}
max=fabs((x[1]-xp)/xp);
for(i=2;i<=n;i++)
if (fabs((x[i]-yp)/yp)>max) max=fabs((x[i]-yp)/yp);
}
while (max>e);
printf("Rozvjazok i perevirka\n");
for (i=1;i<=n;i++)
{
printf(" x[%d]=%2.6f\n",i,x[i]);
printf("f[%d]=%2.6f\n",i,func(x,i));
}
return 0;
}
Результати виконання програми:
x[1]=-0.188679
f[1]=-0.000000
x[2]=0.660377
f[2]=-0.000000
Бло...
